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A Truncated  SQP  Algorithm  for  Large  Scale  Nonlinear 
Programming  Problems  * 

Paul  T.  Boggs  ^ Jon  W.  ToUe  ^ Anthony  J.  Kearsley  § 

August  4,  1992 


Abstract 

We  consider  the  inequality  constrained  nonlinear  programming  problem  and  an 
SQP  algorithm  for  its  solution.  We  are  primarily  concerned  with  two  aspects  of 
the  general  procedure,  namely,  the  approximate  solution  of  the  quadratic  program, 
and  the  need  for  an  appropriate  merit  function.  We  first  describe  an  (iterative) 
interior-point  method  for  the  quadratic  programming  subproblem  that,  no  matter 
when  it  its  terminated,  yields  a descent  direction  for  a suggested  new  merit  function. 

An  algorithm  based  on  ideas  from  trust-region  and  truncated  Newton  methods,  is 
suggested  and  some  of  our  preliminary  numerical  results  are  discussed. 

1.  Introduction 

Large  scale  optimization  problems  are  gradually  submitting  to  the  power  of  ad- 
vanced algorithmic  development  and  of  modern  computing  environments,  leading 
to  the  formulation  of  models  requiring  solutions  of  these  problems  in  a variety  of 
scientific  areas.  Two  excellent  recent  reviews  are  given  by  Coleman  [Col9l]  and 
Conn,  Gould  and  Toint  [ConGT92],  who  survey  some  important  applications  as 
well  as  recent  trends  in  algorithms  and  consider  the  impact  of  parallel  computing 
architectures  for  large  scale  optimization. 

Following  these  authors  we  take  the  term  large  scale  to  mean  any  optimization 
problem  that  is  large  enough  so  that  the  exploitation  of  special  structure  is  im- 
portant. In  this  paper  we  are  particularly  concerned  with  sparsity,  although,  as 
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they  point  out,  other  problem  structures  may  be  important  as  well.  We  assume  the 
general  nonlinear  programming  problem  to  be  of  the  form 

min  f{x) 

X 

subject  to:  g{x)  < 0 

where  / : 3?"  — > 5?^,  and  g : 3?"^.  We  note  that  we  could  include  nonlinear 

equality  constraints  in  {NLP)  without  incurring  any  aiialytical  difficulties,  but  at 
the  expense  of  distracting  technicalities.  We  thus  omit  them  for  the  purposes  of  the 
exposition  here,  but  they  have  been  included  in  our  program. 

Our  basic  tool  for  the  solution  of  {NLP)  is  the  sequential  quadratic  program- 
ming (SQP)  algorithm  in  which,  given  an  approximate  solution-multiplier  pair, 
{x’^,  A*^),  {NLP)  is  approximated  by  a quadratic  program  of  the  form 

min  V/fx*)"*”  -f  B^6 
subject  to:  Vg(x*) ' 6 -f  g < 0. 

Here  5*^  is  taken  to  be  an  approximation  to  the  Hessian  of  the  Lagrangian  for 
{NLP),  i.e.,  for 

^(x.  A)  =1  f{x)  -hff(x)'''A 

we  choose 

In  this  form  the  solution  to  {QP)  provides  a search  direction  for  improving  the 
current  iterate,  x*.  A steplength  is  chosen  in  this  direction  so  cis  to  reduce  a 
merit  function.  Roughly  speaking,  a merit  function  is  a scalar  valued  function 
with  a minimum  at  the  solution  to  {NLP).  Thus  reducing  this  function  ensures 
progress  and  allows  for  the  production  of  a globally  convergent  scheme.  (See  e.g., 
[BogT89]  and  [BogTK9l].)  In  a previous  paper,  [BogTK91],  the  authors  introduced 
a merit  function  for  {NLP)  and  showed  that  it  is  appropriate  for  use  with  the 
SQP  algorithm.  In  this  paper  we  apply  these  ideas  to  the  large  scale  case,  solving 
{QP)  only  approximately  by  an  iterative  interior-point  algorithm  that  we  can  stop 
prematurely.  Such  ideas  are  in  the  spirit  of  truncated  or  inexact  Newton  methods. 
(See  [DemES82]  and,  for  a recent  discussion  of  these  methods,  [EisW91].) 

To  be  more  specific,  we  use  the  interior-point  quadratic  program  solver  of  Boggs, 
et  al.  [BogDRW91].  At  each  iteration  this  method  constructs  a low-dimensional 
subspace  and  solves  {QP)  restricted  to  that  subspace.  We  show  that  halting  this 
procedure  after  any  number  of  steps  yields  a descent  direction  for  the  merit  function. 
The  details  of  this  solver  and  its  properties  relative  to  its  use  in  an  SQP  algorithm 
are  discussed  in  §2. 

The  actual  merit  function  and  a related  approximate  merit  function  are  reviewed 
in  §3.  We  then  state  the  results  just  mentioned,  namely  that  the  inexact  directions 
are  compatible  with  these  functions.  In  §4  we  give  the  details  of  the  algorithm.  One 
of  the  problematic  points  is  how  to  control  the  number  of  iterations  on  {QP)-  Here 
we  use  some  ideas  from  trust  re^ioTi  methods.  We  attempt  to  assess  how  well  {QP) 
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approximates  the  behavior  of  the  merit  function  by  maintaining  an  estimate  of  a 
trust  region  radius.  §4  also  contains  a summary  of  the  results  of  some  numerical 
experimentation  with  the  algorithm  on  a few  large  problems  of  interest.  Our  results 
indicate  that  our  procedure  is  viable  for  large  scale  problems.  Suggestions  for  further 
research  are  contained  in  §5. 


2.  An  Interior-Point  QP  Solver 


Interior-point  methods  for  linear  programming  have  been  demonstrated  to  be  very 
successful,  especially  on  large  problems;  thus  it  is  natural  to  consider  their  exten- 
sion to  quadratic  programs  (QP).  One  method  that  has  performed  well  on  linear 
programs,  and  has  been  extended  to  QP  with  both  good  numerical  results  and 
particularly  interesting  properties  with  respect  to  the  SQP  method  is  the  optimal 
subspace  method  of  Boggs,  et  al.  [BogDRW91].  (See  also  [DomBRW91].)  We  take 
the  QP  of  §1  to  be  of  the  form 


min  6 + 

6 

subject  to:  6 b < 0 


(2.1) 


where  c,  6 € 5?",  Q G A G and  6 G 3?”^. 

The  assumptions  on  (2.1)  that  are  necessary  to  apply  the  algorithm  are  that 
the  problem  is  bounded;  that  A have  full  column  rank;  that  there  exist  feasible 
points  (i.e.,  that  the  constraints  be  consistent);  and  that  Q be  positive  semidefinite. 
Note  that  a full  dimensional  interior  is  not  required.  We  comment  further  on  these 
assumptions  at  the  end  of  this  section. 

Briefly,  the  general  algorithm  can  be  expressed  as  follows. 

OSD  Algorithm  for  Quadratic  Programming 

1.  Given  a feasible  point,  <5°;  set  j :=  0. 

2.  Generate  3 independent  search  directions 

Pit  X = 1, . . . , 3. 

Let  be  the  matrix  whose  columns  are  pi. 

3.  Form  and  solve  the  restricted  quadratic  program 

min  c~^6  -|-  Q6 
subject  to:  A^ 6 4-  6 < 0 

where  6 — 6^  P^Ct  ^Jid  ^ G 3?^.  Call  the  solution  (*. 

4.  Set  :=  6^  -t-  pP^  C*  for  an  appropriate  value  of  the  steplength  p. 

5.  If  stopping  criteria  are  met,  set  J — j,  6j  = 6^  and  exit. 

6.  Go  to  2. 
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The  details  of  the  actual  algorithm  can  be  found  in  [BogDRW91]  and  [Dom- 
BRW91];  here  we  describe  those  that  are  the  most  important  for  its  application 
in  the  SQP  setting.  One  of  the  three  directions  is  always  a descent  direction  with 
respect  to  the  objective  function,  thus  assuring  descent  in  the  objective  value  at 
each  step.  Specifically,  the  algorithm  uses  directions  that  are  solutions  to 


Pi  = U 


(2.2) 


where  ^ is  a positive  scalar  depending  on  the  current  iterate. 


D = diag{l/ri.  A:  = l,...,m}. 


Tk  — —{AkS^  + and  U is  a particular  right  hand  side.  The  form  of  the  matrix 
in  (2.2)  allows  for  efficient  exploitation  of  the  sparsity.  Note  that  if  Q is  positive 
semi-definite,  then  this  matrix  is  positive  definite  for  all  interior  points.  Given  the 
final  direction,  6^,  the  steplength  p is  set  either  to  obtain  the  optimal  solution  in 
the  given  direction  or  to  advance  99%  of  the  distance  to  the  boundary. 

An  important  cispect  of  the  algorithm  is  the  procedure  for  obtaining  an  initial 
feasible  point,  since  we  certainly  do  not  require  that  a feasible  point  be  given.  The 
algorithm  uses  a “Big  M”  method  to  construct  the  Phase  I problem: 


min  c^6  -t- Q 6 + M9 

6,0 

subject  to:  A^ 6 + b — eO  <0 


(2.3) 


where  e is  a vector  of  all  ones  and  6 is  the  “axtificial”  variable.  Clearly  for  6*  large 
enough  the  point  {6,6)  = (0,6*)  is  feasible  for  (2.3).  The  above  procedure  is  thus 
used  until  the  artificial  variable  is  negative,  at  which  point  the  current  value  of  6 
is  feasible,  and  the  M6  and  e6  terms  are  dropped.  If  no  such  value  of  the  artificial 
variable  can  be  found,  then  the  QP  is  not  consistent,  i.e.,  no  feasible  point  exists, 
and  the  algorithm  stops.  In  this  case,  however,  one  can  show  that  the  optimal 
solution  satisfies 

6 = min  max{A7 6 -\-  bj  }, 

6 j ^ 

and  the  resulting  6 is  a reasonable  direction  for  (NLP). 

The  criteria  for  convergence  of  the  algorithm  are  that  at  least  one  of  the  following 
hold:  (a)  the  relative  change  in  two  successive  values  of  the  objective  function  is 
small;  (b)  the  relative  difference  between  the  primal  and  the  dual  objective  function 
values  is  small;  or  (c)  the  relative  difference  between  two  successive  iterates  is  small. 
To  this  list,  we  have  added  the  criterion  (d)  the  scales  solution  vector  exceeds  a 
specified  length.  This  last  condition  has  been  implemented  to  allow  the  use  of  trust 
region  strategies  to  monitor  the  quality  of  the  (QP)  approximation.  In  particular, 
this  procedure  will  cause  the  algorithm  to  halt  if  (QP)  is  unbounded,  again  with  a 
reasonable  direction. 

Note  that  the  assumptions  set  forth  above  ensure  that  a solution  to  (2.3)  exists, 
but  that  the  quadratic  subproblems  arising  in  the  SQP  algorithm  may  not  have 
solutions.  Nevertheless,  the  directions  calculated  by  03D  are  useful  directions  in 
the  solution  of  (NLP). 
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3.  The  Merit  Function 

A merit  function  for  (NLP)  is  typically  a scalar  valued  function  that  has  an  uncon- 
strained minimum  at  x*,  the  solution  to  {NLP).  Thus  a reduction  in  this  function 
implies  that  progress  is  being  made  towards  the  solution. 

In  [BogTK91]  we  derived  a merit  function  for  (NLP)  based  on  the  work  in 
[BogT84]  and  [BogT89]  for  equality  constrained  problems.  This  was  done  by  con- 
sidering the  slack  variable  problem  (see  [Tap80]) 

min  fix) 

1 o (3-1) 

subject  to:  g{x)  -|-  = 0 


where 

S = diag{si,...,s^}. 

The  merit  function  in  [BogT84]  was  then  applied  to  (3.1).  Since  the  resulting  merit 
function  only  contained  references  to  and  not  to  just  Si,  it  was  natural  to  rephrase 
this  merit  function  in  terms  of  Zi  = s?.  This  led  to  the  rather  unusual  situation 
of  having  a constrained  merit  function,  i.e.,  a merit  function  whose  constrained 
minimum  corresponds  to  the  solution  of  (NLP).  Our  merit  function  is 

V’d(x,  z)  = f{x)  + X{x,z)'^c{x,  z)  + ^c(x,  z)A{x,  z)"^c(x,  z)  (3.2) 

where  d is  a scalar, 

c(x,  z)  = g(x)  + Ze 
A(x,z)  = Vg(x)'^Vg(x)  + Z 
A(x,z)  = -A(x,z)~^Vg(x)'^Vf(x) 

and 

Z = diag{zi,...,2,„} 

with  the  z vector  constrained  to  be  nonnegative.  Although  the  merit  function  is 
constrained,  our  algorithm  ensures  that  the  Zi  always  remain  positive;  thus  the 
bounds  present  neither  a theoretical  nor  a computationcd  difficulty.  For  a direction, 
6*^,  in  X obtained  as  the  solution  to  {QP),  we  take  the  direction  for  the  change  in  z 
to  be 

= -[Vff(x*)<5*  + g(x*)-Hz'=]. 

Thus  the  next  step  is 

= X*  -|- 

2*=+!  = z'^  + aq^ 

for  some  value  of  a.  Observe  that  if  2*^  > 0 and  6*  is  feasible,  then 

q*=  + 2*  = -[V3(x*=) V + 3(x*)]  > 0 


(3.3) 
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and  it  follows  that  for  a G (0, 1],  = z*  + aq'^  > 0. 

We  show  in  [BogTK9l]  that  tpd  hcis  certain  desirable  properties  for  d sufficiently 
small.  First,  under  mild  conditions,  a constrained  minimum  of  V’d  corresponds  to 
a solution  of  (NLP).  Furthermore,  the  directions  are  descent  directions 

for  ipd  for  (x*’,z*^)  sufficiently  close  to  feasibility;  a steplength  of  one  is  acceptable 
near  the  solution  if  the  method  is  converging  q-superlinearly;  and  the  directions  are 
always  descent  directions  for  t’(i,z)  ||c|f  (x,z). 

Despite  these  useful  properties,  V’d  has  two  deficiencies  that  preclude  using  it 
directly  in  an  algorithm.  First,  as  stated  above,  6'°  is  only  a descent  direction 
near  feasibihty,  and  second,  it  requires  the  evaluation  of  gradients  and  nontrivial 
computation  to  assess  a prospective  value  of  a.  Thus  we  employ  an  approximate 
merit  function  and  a globalization  strategy  that  overcome  these  deficiencies.  We 
use 

= /(®)  + c(z,2)'^A*=  + ^c(x,z).4*=c(i,z) 

where 


We  show  in  [BogTK91]  that  (S^,  q^)  is  a descent  direction  for  V’d  everywhere,  that 
will  not  interfere  with  rapid  local  convergence,  and  that  the  globalization  strategy 
described  in  §4  is  effective. 

The  main  theoretical  result  described  here  is  that  OSD  and  the  merit  function 
are  compatible.  Specifically,  if  6j  is  only  a partial  solution  to  (QP)  obtained  by 
J iterations  of  OSD  (see  step  5),  the  above  results  continue  to  hold.  We  state  the 
assumptions  that  guarantee  this.  We  use  the  term  strong  local  solution  to  mean  an 
optimal  point,  together  with  a multiplier  vector,  of  {NLP)  at  which  the  following 
hold. 

Al:  The  active  constraint  gradients  are  linearly  independent. 

A2:  Strict  complementary  slackness  holds. 

A3:  The  second  order  sufficient  conditions  hold. 

In  addition  we  make  the  following  assumptions  on  the  {QP)  subproblems: 

A4;  The  matrices  {B*'}  are  uniformly  positive  definite. 

A5:  For  each  k {QP)  at  x*'  has  a strong  local  solution. 

We  also  need  an  assumption  that  guarantees  that  the  merit  function  is  well  defined 
i.e.,  that  A is  nonsingular.  As  in  [BogTK9l]  we  formulate  this  by  partitioning  the 
index  set  of  the  constraints  into  two  subsets  a and  u.  We  can  then  write,  without 
loss  of  generality. 


9{x)  = 


ffa(x) 

gu{x) 
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and  correspondingly, 


Usually  the  index  subset  a will  correspond  to  the  set  of  active  constraints  for  (NLP) 
or  (QP)-  The  necessary  assumption  in  terms  of  a particular  partition  is  the  follow- 
ing. 

A6:  The  set  {Vgt(x)  : i 6 a}  is  linearly  independent  and  > 0. 

A discussion  of  the  implications  of  these  assumptions  for  SQP  algorithms  is  given 
in  [BogTK9l].  The  proofs  of  the  results  make  use  of  the  techniques  in  [BogTK91] 
combined  with  an  induction  argument. 


4.  Algorithm  and  Num.erical  Results 

A brief  statement  of  the  final  algorithm  is  as  follows.  Following  the  statement,  we 
give  a brief  discussion  of  some  important  points. 

SQP  Algorithm 

1.  Given  x°,  r (trust  region  radius),  t]  (globahzation  parameter),  and  d (merit 
function  parameter): 

Set  k :=  0. 

2.  Using  OSD,  iterate  while  ||6||  < r on 

minV/(x*)'''^  + 

6 

subject  to:  V3(x*')^<5  -|-  g{x^)  < 0 

to  obtain  6*. 

3.  Set  = -[Vg{x’^)^6^  + g{x^)  + 2*]. 

4.  (Globalization  Step) 

Choose  a*  such  that  is  reduced. 

If  ||c(x*  -I-  > ||c(5c*', -z*^)||  and  ||c(x*, 2*)||  > t], 

reduce  q*  until  ||c(x* -|- 2* -f  a*=9*)||  < ||c(x^,  2*)|j. 

5.  If  V’(i(a:*’  -f  a^6^,z^  -|-  a^q^)  > V’<i(a:*',  2*) 
set  T]  = ^ ||c(i*^, 2*')||. 

6.  Set 


:=  -f 

2*+i  :=  2*+Q*g*. 


7.  If  convergence  criteria  are  met,  quit. 

8.  Adjust  T. 
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9.  Set  k k + 1]  goto  2. 

A few  comments  are  necessary.  First,  the  globalization  step  is  based  on  the 
work  in  [BogT89].  In  brief,  77  is  an  estimate  of  the  radius  of  the  domain  containing 
the  feasible  region  in  which  the  true  merit  function,  is  reduced  in  the 

direction  {6,q).  For  all  iterates,  the  algorithm  first  requires  that  the  approximate 
merit  function  be  reduced.  If  the  current  iterate  lies  outside  the  77-domain,  then  the 
algorithm  also  requires  that  the  constraint  infeasibilities  be  reduced.  If  the  iterate 
lies  inside  the  77-domain  then  the  true  merit  function  should  also  be  reduced;  if 
not,  then  77  is  reduced.  This  allows  steps  that  may  increase  the  merit  function,  but 
only  in  a controlled  way.  The  steps  that  increase  the  merit  function  are  usually 
seen  only  in  ecU'ly  iterations  or  after  active  set  changes.  Second,  our  procedure  for 
updating  r is  to  compute  the  predicted  relative  reduction  of  the  merit  function 
based  on  the  (QP)  and  compare  that  with  the  actual  relative  reduction.  This 
comparison  of  predicted  and  actual  reductions  is  done  using  the  approximate  merit 
function  if  the  current  iterate  lies  outside  the  77-domain.  The  true  merit  function 
is  employed  otherwise.  We  then  use  standard  updating  strategies  to  adjust  r (see 
e.g.  [DenS83]  or  [MorS83]).  Third,  the  penalty  parameter,  d,  is  updated  in  a very 
straight  forward  manner.  Essentially,  an  estimate  of  the  condition  of  the  problem 
is  monitored.  In  the  event  that  this  estimate  increases  significantly,  d is  decreased. 
Provided  that  the  initial  value  of  d is  reasonable,  this  updating  did  not  occur  often, 
and  is  only  observed  when  the  iterates  are  outside  the  77-domain.  The  algorithm 
did  not  ‘hug’,  or  stick  too  closely  to  the  constraint  manifold,  as  is  the  case  when 
the  penalty  parameter  becomes  too  small.  Computationally,  this  simple  procedure 
for  updating  d appears  to  be  effective  even  in  the  presence  of  highly  nonlinear 
constraints. 

We  have  used  this  procedure  to  solve  several  problems  in  the  range  of  100-500 
variables  with  up  to  500  constraints.  Many  of  these  problems  have  arisen  from 
discretizations  of  control  problems  where  the  Hessian  of  the  Lagrangian  and  the 
Jacobian  of  the  constraints  have  some  known  sparsity  structure.  These  problems  are 
somewhat  special,  in  that  we  knew  that  the  major  expense  in  the  calculation  of  an 
iterate  comes  from  the  solving  of  the  {QP)  . Typically,  the  constraints  are  nonhnear 
inequalities  that,  in  some  way,  hmit  the  control  variables,  and  the  objective  function 
is  an  energy  approximation.  The  number  of  constraints  is  greater  than  the  number 
of  variables  in  many  of  the  problems  we  solved.  In  our  testing,  we  use  forward  finite- 
difference  approximations  to  gradients  and  Hessians,  and  modify  the  Hessian  of  the 
Lagrangian  to  be  positive  semidefinite  in  cases  where  it  is  not  (e.g.  [GilMW81]). 
This  latter  procedure  requires  the  addition  of  a non-negative  diagonal  matrix  to  the 
Hessian  approximation. 

Our  observations  include  the  following. 

• The  number  of  major  iterations  is  reasonable. 

• The  globalization  procedure  remains  efficient,  i.e.  many  full  steps  are  ac- 
cepted. 

• Close  to  the  solution,  the  trust  region  becomes  inactive. 
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• The  trust-region  strategy  is  basically  effective,  i.e.,  it  prevents  long,  unprof- 
itable steps  from  being  generated  at  the  beginning  and  after  the  active  set 
changes. 

• In  general  a small  number  of  iterations  of  OSD  suffices  at  each  major  iteration, 
and  a very  small  number  near  the  solution. 

• As  in  all  of  our  previous  work  in  this  area,  the  parameter  d in  the  merit 
function  is  not  critical,  i.e.,  the  performance  of  the  algorithm  is  not  changed 
much  by  changes  in  the  strategy  for  adjusting  d. 

5.  Future  Research 

We  have  described  a preliminary  version  of  an  extension  of  the  SQP  ideas  to  the 
large  scale  case.  In  doing  so,  we  have  used  a combination  of  an  interior-point 
method  for  solving  (QP)  with  trust  region  and  truncated  Newton  methods  to  create 
a promising  algorithm.  There  are,  however,  many  computational  and  theoretical 
aspects  of  this  algorithm  that  need  further  analysis  and  testing.  Computationally, 
we  need  to  continue  to  test  the  procedure  to  discover  its  strengths  and  weaknesses. 
At  the  same  time,  the  limitation  on  the  QP  solver,  OSD,  that  the  Hessian  must  be 
positive  semidefinite,  appears  to  be  surmountable.  In  particular,  we  believe  that 
the  solution  of  the  reduced  quadratic  program  (step  3 of  OSD)  can  be  modified  to 
handle  an  indefinite  (or  negative  definite)  Hessian.  This  would  allow  us  to  avoid  the 
extra  work  of  ensuring  that  the  Hessian  is  positive  definite,  and  to  explore  directions 
of  negative  curvature. 

Our  theoretical  analysis  described  in  §3  relies  on  the  usual  strong  assumptions 
that  are  typically  satisfied  in  the  small  scale  case.  Some  of  these  assumptions, 
however,  are  often  not  satisfied  in  large  problems.  In  particular,  large  problems  may 
be  highly  degenerate.  We  know  that  the  interior-point  algorithms  for  LP  and  QP 
have  no  problem  with  these  cases.  Some  of  the  problems  that  we  have  attempted 
have  been  degenerate,  and,  although  the  theory  does  not  apply,  the  algorithm 
had  no  difficulty  in  solving  them.  Also,  in  some  of  the  problems,  the  queidratic 
subproblems  were  not  always  consistent.  This,  too,  caused  no  difficulty  for  the 
algorithm,  but  is  a problem  for  the  theory.  Thus,  obtaining  good  theoretical  results 
under  a weakened  set  of  assumptions  is  an  important  task  for  further  research. 
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